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This paper addresses two problems of interest in service 
system analysis: (a) that of making statistical, data- 

driven estimates of the long-run probability of a long 
delay, and (b) the assessment of rate of approach to a 
long-run system performance measure such as expected 
delay, the rate being characterized by a simple exponent 
tial , at least initially. Both are illustrated by 
reference to M/G/l and related systems. 



1 . Introduction 

i 

The application of probability theory to a wide variety of* congestion problems 
arising in communication systems has been well catalogued in many papers; the 
treatises by Borovkov (1984), Cohen (1969), Cooper (1972), Cox and Smith (1961), 
Feller (1966), Kleinrock (1975), Gross and Harris (1974), and many others attest 
to the attractiveness of the area for probability modelling. In this body of work, 
several features are noticeable. First, most of the elegant solutions obtained 
are in somewhat implicit form, being presented as functional equations, or, fre- 
quently, as integral (Laplace) transforms, generating functions, and sometimes as 
combinations of the above. Moments (particularly the first) of delays or number 
in queue under long-run conditions when steady state prevails are the most explicit 
and assessible performance summaries. Second, the results obtained are presented 
in terms of component distribution functions and stochastic processes (renewal, 
Poisson, etc.) that are taken as known; only rarely are issues addressed that 
arise when actual data is to be used as a basis for inference from the models; 
however, see Cox (1965). Third, only fragmentary comprehensible information con- 
cerning transient behavior, not to mention time-dependence of parameters is avail- 
able; recently however work by Roth and Odoni (1983), Roth (1981), Abate and 
Whitt (1985), Lee (1985) has elucidated the simple exponential -1 ike approach to 
the steady state displayed by a great many stochastic systems with stationary, 
time-homogeneous arrivals and service processes. 

Instead of directly dealing with the above classes of problems much attention 
has been concentrated on modelling large networks of servers, particularly steady- 
state Markovian, cf. Kelly (1979) for an elegant treatment; see package programs 
created by AT&T (the D. Mitra PANACEA) and by IBM (work by Lavenberg, et al ) . 

Other numerical work by Neuts (1981) and associates has pioneered the exploita- 
tion of special structure of certain stochastic systems. And a considerable ef- 
fort to create approximations, ranging from heavy-traffic and diffusion to, lately, 
WKB approaches, see Knessl , Matkowsky et^ aj_ (1986) has been evident; see in par- 
ticular Newel-1 (1982) for pioneering work. 

This paper deals specifically with approaches to the two classes of somewhat 
neglected problems referred to earlier: those of data-driven inference, and of 

assessing transient behavior. The approaches proposed are illustrated concretely 
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in terms of the simple M/G/l system, but apply more widely. 

2. Inference Concerning Long Steady-State Delays in M/G/l -Like Systems 



Consider a single service system approached by stationary Poisson (A) traffic 
with A known. Service times, or message lengths, X, are independently distributed 
as Fx(x); assume A E [X] = p < 1. This is the M/G/l system. Let observations of 
the service times be all that is known about F: denote these by x-j ,X 2 » • • • »x n . 

The objective is to supply estimates of long-run characteristics of the system: 
in particular, estimates, and error assessments thereof, of the probability of a 
long delay experienced by an arriving message. 

2.1 Virtual Waiting Time Transform; The Long Tail. 

It is well known that if W(t) is the virtual waiting time in the M/G/l and 
p < 1, then the Laplace transform 

E[e" sW ] = lim E[e" sW(t ^] = 1 ~ p ' ■ ■ v = — - £ — (2.1) 

t-~ Thrill [el'll ] 1-PA(S) ** - 1 

1 PL sE[x] j 

where A(s) is the Laplace-Stiel tjes transform of a distribution. It can also be 
seen that in various interrupted and priority service systems, including those in 
which the server "takes vacations," that the above can be changed to 

E[e" sW ] = ] : p • C(s) (2.2) 

l-pB(s) 

where possibly B(s) is the transform of a completion time (effective service time, 
i.e. X including durations of a random number of interruptions occurring therein), 
and C(s) is the transform of an honest distribution that accounts for effects of 
busy period starts; see Gaver (1968). 

If A(s) (or B(s)) exists for s > s 0 , s 0 < 0 then there will be a smallest 
real zero, -s = < > 0 , of the denominator of (2 . 1 ),( ( 2 . 2 ) ) , which can be used to 
show that 



P{W > w } ~ DMe"™, w oo (2.3) 

a 

One way of establishing the above exponential tail property is to introduce 



ip (s) 



sW- 



. L -Ete-n a 



/ P{W >w}e~ SW dw 



into ( 2 . 1 ) which leads to 

>Ks) = p[- " 3 ^ ] + pA(s)i|/(s) 

equivalent to the terminating renewal equation, see Feller (1966), p. 362, 



(2.4) 



Fi» ( w ) = P{W >w} = H(w) + p / F(w-x)A(dx) 

0 

Introduce F^(w) = F w (w)e 0W , 9 real and pos-itive, into (2.5) 

Fu(w) = H # (w) + / F^(w-x)pe 0 X A(dx) , 

0 



and choose 9 = < so that 

oo 

9x 



/ pe 0 X A(dx) = / A # (dx) = 1 , 

0 0 



(2.5) 



( 2 . 6 ) 



(2.7) 
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—4 

yielding a standard renewal equation for F , to which the key renewal theorem 
applies, yielding as w °° 



lim F^(w) = 1 im F^(w)e kW 

W*^°° 



D(K) 



/ H^(w)dw 
0 

00 it 

/ x A (dx) 

0 



( 2 . 8 ) 



equivalent to (2.3). A similar expression is available for (2.2). 

2.3 Statistical Estimation of Probability of Long Wait. 

If data are available on service times (message lengths) then an estimate of 
the requisite transform is 



~ „ v i n -sx. 

E[e' SX ] • £ I e ' 

1=1 

and k is the unique positive solution of this sample equivalent of (2.7): 



n ex . 

I 

i=l 



r "- 15 =1 



(2.9) 



( 2 . 10 ) 



A three-term Taylors' series expansion around 0=0 gives an initial estimate 

k = 2(— ) (2.11) 

7? 

nr k k 

where x = (x-j + ... + * n )/n, the kth_ sample moment; (2.11) is recognized as being 

the non-parametric sample version of the exponential distribution parameter ob- 
tained by normalizing to W/E[W] and letting p t 1 (the heavy traffic approximation). 
A search or Newton-Raphson iteration starting from (2.11) quickly yields k as 
accurately as is necessary. 



Some theoretical asymptotic results concerning modes of behavior of < will now 
be sketched; more detail will be presented elsewhere. First, it Is essential that 

E[e K ^] < °° in order for (2.7) to hold, and so by continuity f^ m ^(<) = E[X m e K ^] <°° 
for any integer m. Now express (2.10) as 

~ i 1 n 0 x . 

f(«) - l e 1 -1)}- 1 (2.12) 

6 n 1=1 

/s /s 

and expand in Taylor's series around <: since f(<) = 0, for some $ 

0 = f(<) + (K-K)f'(ic) + |-(k - <)^f " ( qx) , (2.13) 



so 



K - K 



f(«)^ 

f ' (k) + ^( k-k) f " ( (J ik) 



'lote that if E[e A ] < <» then Var[e A ] < °° and 

2 

Var[ft<)3 = Vjf Var[e KX ] , 

< 

vhi 1 e 



(2,14) 



(2.15) 
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E[f(<)] = 0 . 



Additionally, f'(K) -*■ E[f’(i<:)] by the weak law of large numbers and <-< -*■ 0 in 
probability, so by the central limit theorem. 



of the sum f(<) must be in the domain of attraction of a stable law, which will be 
assumed. In any case ?^(k) -*• E[?'(k)] as before, so a proper normalization by a 
power of n shows that (k-k) has a stable law as limiting distribution. Fortunate- 
ly, the relatively well-behaved Normal/Gauss limit prevails when traffic intensity 
is relatively high (p > 1/2 in case X ~ Exp(p), for example). - _j 

2.4 Assessing Variability of the Estimates. 

Having obtained an estimate of the parameter k, and of the probability of a 
waiting time exceeding large t, which involves k, it is desirable to appraise the 
errors involved. Two error sources are: the systematic error -(bias) resulting 

from fitting an incorrect model, and the effect of finite sample size (random 
error). The latter is the easiest to evaluate, provided the underlying model is 
nearly correct. The bias issue is more difficult to deal with; one approach is 
to begin by fitting an elaborate, mul ti parameter model and then to check for the 
contribution of extra parameters in a prediction context. Since this paper takes 
a non-parametric approach to estimation, the bias resulting from an incorrect 
specification of the service time distribution is presumably absent, although the 
assumptions of stationary Poisson arrivals, and independence are still reflected 
in the basic transform (2.1), and hence in the estimates. Unfortunately, classi- 
cal procedures that, for example, avail themselves of the asymptotic approxima- 
tions of properties of maximum likelihood parameter estimates to estimate sampling 
errors are no longer available in the present environment. We have chosen in- 
stead to investigate the performance of several modern procedures often recom- 
mended for obtaining standard errors of estimate and confidence Intervals: the 

jackkni fe (Quenouille, Tukey, Miller, Hinkley, and others) and the bootstrap 
(Efron and associates and others). These methods, particularly the bootstrap, 
which involves simulation, are computer intensive, but are the only alternatives 
currently known for dealing with complex situations such as that at hand. 

The Jackknife 

■ l 

The above name refers to a procedure originally introduced by Quenouille for 
bias reduction (1956), and adapted by Tukey (1958) to obtain approximate confidence 
intervals. Suppose interest is. on a parameter 6 (e.g. our k, or some function 
thereof) that is estimated by 6, using a complex calculation from data 
(xi , x 2, x . . »x n ) » just as k is. The idea is that of assessing variability by recom- 
puting k after removing independent subgroups of data of equal size, and then 
using the recomputed k values to estimate a variance, which is in turn applied to 
state a standard error or a two-sided confidence interval that contains the true 
9 with specified confidence. A few details follow; for more, see Efron (1982) 
and his more recent work, or Mosteller and Tukey (1977). The actual calculation 
involves splitting n into g disjoint groups of size m; n = mg. Then calculate 
9(-j)> j = l,2,...,g: the estimate of 0 that omits the jth[ group. Now Tukey 

(but not Efron) computes pseudo-val ues 



/n(<-<)* — ^ ^ ~ N(0,1) 




(2.16) 



2kx 

the standard Normal /Gaussian distribution. If, however, E[e ] = °°, then if any 
kind of limiting distribution is to exist, the components 



yj = ge - (g-i)e ( _ J) 
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which are then treated as independent: use y = 0 JK or the point estimate of 0, 

and its approximate variance as 

Sy/9 - j ( y j -y ) 2 / ( g - 1 ) ( g ) * j < Vj)' ? (-j ) )2(9 ' 1 )/9 ' 

J 1 j 

as it turns out. Tukey recommends referring y to Student's t with g-1 degrees of 
freedom to obtain confidence limits. As the name "jackknife" is intended to 
imply, the tool is inexact and a bit crude for small samples--just as a true 
jackknife is not well-adapted for delicate surgery. It is a handy non-parametric 
option . 

The Bootstrap 

In simplest form the bootstrap suggests creating from observations (x., 
i = l,...,n) the empirical distribution 

i 1 n ( 1 , x. £ x 

F n (x) ° n X. 1 < x i ,x); '<V X > ’ j 

i=l ( 0, x. > x 

i 

Then one re-samples wi£h replacement from F n to obtain a bootstrap sample ofsize 
n, from which 0 (e.g. k( 1)) is calculated. This^process is repeated B (e.g. 
300-500) times, and the empirical distribution, F~, of the resulting estimates, 
{<(b), b = 1,2,..,B) is employed as an approximation to the sampling distribution 
of k: the upper and lower 10% points of F~ approximate two-sided 80% confidence 
limits for k, for instance. K 



3. Transient Response: Exponential Characterizations 

Applications of queueing theory frequently require assessment of service sys- 
tem transient behavior, particularly when the system is initially idle and heavy 
traffic conditions prevail, for then approach to steady-state conditions may be 
slow. Odoni and Roth (1983) have reviewed work on the nature of the approach to 
steady-state conditions, and have investigated characterizations of that approach 
that are simply stated and comprehensible, being exponential . In this section 
we elaborate upon the theme of exponential approach. Work by Abate and Whitt 
(1985) has similar aims, but proceeds somewhat differently. 



3.1 From Transform to Exponential 

Suppose (X(t), t >0} is the state variable of a process, and that expressions 
for 



oo oo 

(a) Us) = j e' St E|>(X(t)) j X(0) ]dt = / e" s V(t)dt 

A 0 0 A 

(b) 1 im E[ip(X(t) | X(0) : < °° 

are known. An approximate representation of 

(c) <p x (t) = E[>(X(t)) j X( 0)] 



(3.1) 



Is desired. The proposed representation is a simple exponential interpolation of 
the form 



if» x (t) = ^ x (0)e" Bt + <Jj x H( 1 -e' Bt ) , 



(3.2) 



where (3 governs the rate of approach to the steady-state value. The proposal is 

5 



(3.3) 



to identify 3 by minimizing the integrated squared error 

oo 

MB) = / 4 x (t) -^ x (0)e" 3t +^ x (~)(l -e" 6t )} 2 w(t)dt . 

The weight function w(t) may be made to focus upon values of g appropriate for 
t-regions of interest. To optimize on B, study 

oo 

^dB^ = K J { ^X (t) -^ X (0)e ’ 6t ' ^x (o ° )(1 -e" 6t )}te" 3t w(t)dt (3.4) 

= 0 . 

This is equivalent to 

oo oo oo 

/ ip Y (t)w(t)te ^dt = ip Y (0) / e ^tw(t)dt + iJj y («>) / (1-e ^)e ^ t tw(t)dt (3.5) 
0 0 0 

If w(t) = 1 the above can be expressed in terms of the derivative of jthe Laplace 
transform ^(s): 



= ~T^X^ + T ’ (3.6) 

6 

often a transcendental equation that must be solved numerically for 3. Clearly, 
if no solution or multiple solutions exist then the simple universal form (3.2) 
is called into question, but in such cases weight functions are usefully invoked. 
If interest centers on ultimate approach (t -*• °°) to ^„(°°) then the smallest 3 
satisfying (3.6) is of interest. 

3.2 The M/G/l Example. 

As an important illustration of the above ideas, let X(t) = W(t) be the 
virtual waiting time in an M/G/l system, and consider (J^j(t) = E[W(t)|W(0) =0]. 

We know that 



where 



M s) = + P 00 (s) T 



p 00 (s) = s+A[l-b(s)] 



UX [kMiI] 

s 



ih 

s 



(3.7) 



(3.8) 



b(s) being the transform of a busy period duration, satisfying 

b(s) = F x (s + A( 1 -b( s) ) ; 

/V 

Fx is the Lapl ace-Stiel tjes transform of the service time or message length dis- 
tribution Fx-' Numerical solution of (3.6) is now possible. An approximation to 
the (smallest) value of 3 governing the ultimate approach (t -*■«>) when p t 1, i.e. 
in the heavy-traffic situation, is 



PECX 2 ] 

1 -P 



{- 



1 



•h"’ (0) 



-} 



(3.9) 



In general, h"'(0) will involve higher moments of the service time; in case 
X ~ Exp(p) considerable simplification occurs and 



3 - p(l-p) 3 {8(1 +p)} 3,3 
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(3.10) 



vhich will not, of course agree exactly with formulas or numerical results given 
previously by other authors, e.g. see Odoni and Roth (1983) being a compromise 
["Procrustean" ) interpolation of a single exponential over the infinite t-range. 
[n order to focus upon a 8-value relevant to ultimate approach it has been found 
:hat a weight procedure of the form w. + ,(t) = [l-exp(-8(k)t)] is workable and 
tractable; start with 8(0) = 0 (the unweighted procedure described), then intro- 
Juce w^ to determine 8^ and iterate to convergence, which occurs rapidly. 

\ similar approach should apply to find a 8 appropriate for small t. 



The procedures described above should be applicable in the data-driven context 
)f the problem of Section 1: a first step is to introduce the empirical transform 

and from it proceed to b(s), to tj) w (s), and then to 8* Progress in this area, 
ind details of the above investigations will be reported elsewhere. 
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